################################################################
### Generate illustrative examples for "Long-Term Effects for Models with Temporal Dependence"
###
### Original code by dave carter & curt signorino, 4/27/2007
################################################################



############################### Negative (weibull) ###############################################
        Nobs<<-2000
        maxperiods=140
        t=1:(maxperiods+10)
        
	# constant and b_x
	  b0 = -3
	  b1 = 1

        # decreasing failure rate
	   scalx=2
	   locx=.5
	   lam = 1
         h = lam*(.25)*(lam*t)^(.25-1)
         h = h[2:length(h)]*12

       # generate data

        x=-999
        y=-999
        tim=-999
        group=-999
        gg=1
        tt=1
        for (i in 1:Nobs){
          xx=(runif(1)-locx)*scalx         # uniform dist covariate
          #p=1-exp(-exp(-3+xx+h[tt]))       # cloglog link 
          p=1/(1+exp(-(-3+xx+h[tt])))      # logit link
          yy=(runif(1)<p)                  # y is binary
          x=rbind(x,xx)                    # append x
          y=rbind(y,yy)                    # append y
          tim=rbind(tim,tt)                # append time
          group=rbind(group,gg)            # append group
          tt=(yy==0)*(tt+1)+(yy==1)*(1)    # update time counter
          gg=(yy==0)*(gg)+(yy==1)*(gg+1)   # update group counter
        }

        x=matrix(x[2:(Nobs+1)],Nobs,1)
        y=matrix(y[2:(Nobs+1)],Nobs,1)
        tim=matrix(tim[2:(Nobs+1)],Nobs,1)
        group=matrix(group[2:(Nobs+1)],Nobs,1)

	gg<-as.numeric(x)
	df = data.frame(x, y, tim, group)        

	write.dta(df, "NegativeFF.dta")

############################### Positive (weibull) ###############################################
        Nobs<<-2000
        locx=.5
        scalx=2
        maxperiods=50
        lam=0.3
        t=1:(maxperiods+10)

        # increasing failure rate
         h = (lam/5*t)^(1.5)

       # generate data
        x=-999
        y=-999
        tim=-999
        group=-999
        gg=1
        tt=1
        for (i in 1:Nobs){
          xx=(runif(1)-locx)*scalx         # uniform dist covariate
          #p=1-exp(-exp(-3+xx+h[tt]))       # cloglog link 
          p=1/(1+exp(-(-3+xx+h[tt])))      # logit link
          yy=(runif(1)<p)                  # y is binary
          x=rbind(x,xx)                    # append x
          y=rbind(y,yy)                    # append y
          tim=rbind(tim,tt)                # append time
          group=rbind(group,gg)            # append group
          tt=(yy==0)*(tt+1)+(yy==1)*(1)    # update time counter
          gg=(yy==0)*(gg)+(yy==1)*(gg+1)   # update group counter
        }

        x=matrix(x[2:(Nobs+1)],Nobs,1)
        y=matrix(y[2:(Nobs+1)],Nobs,1)
        tim=matrix(tim[2:(Nobs+1)],Nobs,1)
        group=matrix(group[2:(Nobs+1)],Nobs,1)

	gg<-as.numeric(x)
	df = data.frame(x, y, tim, group)        

	write.dta(df, "PositiveFF.dta")

############################### Non-Monotonic (parabolic) ############################################
        Nobs<<-2000
        locx=.5
        scalx=4
        maxperiods=100
        lam=0.25
        t=1:(maxperiods+10)

        # non-monotonic parabolic failure rate
         h = 1/3*(lam/2*(t-16))^2
	   h = h*2 #multiplying hazard to mitigate problem of scale relative to x's

       # generate data
        x=-999
        y=-999
        tim=-999
        group=-999
        gg=1
        tt=1
        for (i in 1:Nobs){
          xx=(runif(1)-locx)*scalx         # uniform dist covariate
          #p=1-exp(-exp(-3+xx+h[tt]))       # cloglog link 
          p=1/(1+exp(-(-3+xx+h[tt])))      # logit link
          yy=(runif(1)<p)                  # y is binary
          x=rbind(x,xx)                    # append x
          y=rbind(y,yy)                    # append y
          tim=rbind(tim,tt)                # append time
          group=rbind(group,gg)            # append group
          tt=(yy==0)*(tt+1)+(yy==1)*(1)    # update time counter
          gg=(yy==0)*(gg)+(yy==1)*(gg+1)   # update group counter
        }

        x=matrix(x[2:(Nobs+1)],Nobs,1)
        y=matrix(y[2:(Nobs+1)],Nobs,1)
        tim=matrix(tim[2:(Nobs+1)],Nobs,1)
        group=matrix(group[2:(Nobs+1)],Nobs,1)

	gg<-as.numeric(x)
	df = data.frame(x, y, tim, group)        

	write.dta(df, "Non-Monotonic1FF.dta")



############################### Non-Monotonic (log-logistic) ############################################
        Nobs<<-2000
        locx=.5
        scalx=4
        maxperiods=600
        lam=0.125
        t=1:(maxperiods+10)

        # non-monotonic failure rate
         h = (lam*(1.5)*(lam*t)^(1.5-1))/(1+(lam*t)^(1.5))
	   h = h*30 #multiplying hazard to mitigate problem of scale relative to x's

       # generate data
        x=-999
        y=-999
        tim=-999
        group=-999
        gg=1
        tt=1
        for (i in 1:Nobs){
          xx=(runif(1)-locx)*scalx         # uniform dist covariate
          #p=1-exp(-exp(-3+xx+h[tt]))       # cloglog link 
          p=1/(1+exp(-(-4+xx+h[tt])))      # logit link
          yy=(runif(1)<p)                  # y is binary
          x=rbind(x,xx)                    # append x
          y=rbind(y,yy)                    # append y
          tim=rbind(tim,tt)                # append time
          group=rbind(group,gg)            # append group
          tt=(yy==0)*(tt+1)+(yy==1)*(1)    # update time counter
          gg=(yy==0)*(gg)+(yy==1)*(gg+1)   # update group counter
        }

        x=matrix(x[2:(Nobs+1)],Nobs,1)
        y=matrix(y[2:(Nobs+1)],Nobs,1)
        tim=matrix(tim[2:(Nobs+1)],Nobs,1)
        group=matrix(group[2:(Nobs+1)],Nobs,1)

	gg<-as.numeric(x)
	df = data.frame(x, y, tim, group)        

	write.dta(df, "Non-Monotonic2FF.dta")

